Method and system for simplifying models

ABSTRACT

An example method of simplifying a system expressed as differential algebraic equations includes transforming the differential algebraic equations into Hessenberg form, the Hessenberg form including algebraic equations and differential equations, using the algebraic equations to express a constraint manifold, and using a combination of differentiation and projections onto the normal and tangential spaces of the constraint manifold to simplify the differential equations.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority to U.S. Provisional Application No.61/415,094, filed on Nov. 18, 2010.

FIELD OF THE INVENTION

The present invention relates to a method for simplifying a modelexpressed as a system of differential algebraic equations (DAE).

DESCRIPTION OF THE PRIOR ART

Systems may be modeled as mathematical expression to enable the systemto be tested. Modeling systems is particularly useful during the designprocess as these mathematical expressions allow the system to be testedbefore a prototype is produced. Specifically, modeling the system allowsan engineer to see how a system performs, for example mechanical systemssuch as engines under various loads and under environmental factors.

Systems may be expressed symbolically using variables which show how thesystem behaves. This allows an engineer to test how the system respondsunder certain conditions. Thus, the system may be subject to variousconditions before a prototype is built. The conditions may be expressedin variables which represent the behavior of the system while othervariables represent forces on the system which enforce the conditions onthe system.

For example, some of the variables may represent external forces such asfriction, heat, and the like. Accordingly, these variables may depend onother factors such as time, or another condition of the model. Thus,taking into account the internal behavior of the system and the externalforces acting on the system, the system may be expressed in a system ofdifferential algebraic equations.

The expression of a system may include numerous algebraic variables anddifferential variables. Processing such an expression may be timeconsuming. Accordingly, it is desirable to reduce the number ofalgebraic variables and differential variables within the system,simplifying the system.

Accordingly, it is highly desirable to have a method of reducing thenumber of algebraic and differential variables; and algebraic anddifferential equations, in a mathematical model of a system so as toreduce the processing time required to solve these systems numerically.

SUMMARY OF THE INVENTION

Examples of the present invention relate to methods for simplifyingmodels expressed as a system of differential algebraic equations (DAEs).Reduction of model represented by DAEs to a system of ordinarydifferential equations (ODEs) facilitates solving the model numerically.Once ODEs are obtained, any appropriate ODE solution approach can beused to obtain a numerical solution of the model. Examples of thepresent invention simplify the reduction of DAEs to an ODE system, e.g.by reducing the number of differential and/or algebraic variables.

Examples of the present invention include combining index reduction withprojection onto the constraint manifold. Index reduction benefits fromthe projection approach for DAEs of various forms, in particular DAEsthat appear in the application of optimization under constraints.Projecting the model onto the normal and tangent spaces of theconstraint manifold allows the elimination of algebraic and/ordifferential variables, simplifying the model.

Applications include simplification and numerical solution of models ofvarious physical systems, such as models of mechanical systems such asengine models, models of electrical systems such as electrical circuitmodels, and the like. For any physical system, mathematical models arederived based on the relevant physical laws. Constraints may representnatural conditions imposed on the system by the physical laws.

In a computational environment including at least one processor, anexample method of simplifying a model of a physical system expressed asDAEs includes differentiating the algebraic equations to produce basesfor a normal space and a tangent space of the constraint manifold, andreducing the number of differential variables and/or algebraic variablesin the model using index reduction of the differential equationscombined with projection of the differential equations onto the normaland tangent spaces.

In engineering development, prototyping and testing can be replaced byvirtual modeling and simulation. A physical system can be represented bya mathematical model within a software modeling and simulationenvironment, and such models often include DAEs that are solvednumerically. Simulation can be achieved on a computer for variousscenarios. Improved simulation approaches allow more rapid developmentand prototyping. In the early stages at least, prototyping usingphysical objects (such as engine components or electronic components)can be replaced by virtual modeling. Hence, improved simulations allowmore rapid development of products or processes. Any approach thatreduces the computational demands of finding a numerical solution, andhence reduces simulation times, allows more extensive use of virtualmodeling. This reduces product development times, and hence improvednumerical solution approaches are of considerable value and interest ina variety of technical fields.

BRIEF DESCRIPTION OF THE DRAWINGS

FIGS. 1A-1B illustrate an example method of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

Examples of the present invention relate to a method for simplifying amathematical model of a system. Computer simulations of physical systemsoften require the numerical solution of a model. Simplification of DAEsto ordinary differential equations (ODEs), which can then be solvednumerically, is very useful for dynamical modeling and simulationapplications. Applications include improved modeling of physicalsystems, such as constrained mechanical systems, electrical circuits,models generated using high level modeling tools, and optimizationproblems constrained by one or more physical laws. However, theapproaches described are not limited to any particular system.

An example method for numerically solving a model in the form of adifferential algebraic equations (DAEs) includes expressing the DAEs inHessenberg form, providing representations of the algebraic equations asa constrained manifold, determining the normal and tangential spaces ofthe manifold, and projecting the system onto the tangential and normalspaces of the manifold so as to reduce the number of algebraic anddifferential variables in the expressed system.

A DAE system can be simplified by projecting the system onto normal andtangential spaces of the constraint manifold, and solving the projectedsystem symbolically so as to reduce the number of differential andalgebraic variables within the DAEs.

According to an example of the present invention, a method of modelsimplification include transforming the DAE into Hessenberg form,considering the algebraic equations as a constraint manifold, findingthe normal and tangential spaces to the constraint manifold, projectingthe system onto the normal and tangential spaces, where possible solvingalgebraic equations found by the normal projection, and using thetangential projection to give ODEs with fewer unknowns than the originalsystem.

A DAE expressed in its most general form may not be immediately suitablefor simplification using a projection method. However, expression inHessenberg form allows projection and index reduction to be combined.

Examples of the invention may be achieved using a computer, including aprocessing unit and associated components, for example executing aHessenberg transformation algorithm operable to transform the DAE intoHessenberg form.

Hessenberg Index-1

The general form of a Hessenberg Index-1 may be expressed as follows:{dot over (x)}=f(t,x,z)0=h(t,x,z)  (1)wherein x=x(t)=(x₁(t), . . . , x_(n)(t))^(T) is the vector ofdifferential variables f=(f₁, . . . , f_(n))^(T), and z=z(t)=(z₁(t), . .. , z_(k)(t))^(T) is the vectors of algebraic variables h=(h₁, . . . ,h_(k))^(T), and the Jacobian ∂h/∂z is nonsingular for all t.

In what follows, we consider a Hessenberg form, which is reduced to anindex 1 by differentiation. The reduced Hessenberg form, meaning hereina Hessenberg form reduced to an index of 1, as in Equation (1), is thenexamined to determine if projecting the system will result in areduction of algebraic variables and differential variables. Projectingor projection used herein describes the process whereby the system isprojected onto the normal and tangential spaces of the constraintmanifold. Naturally, some systems may have a higher index than others.Accordingly, each system must be evaluated to determine if projectionwill reduce the number of differential and algebraic variables of thesystem.

Hessenberg Index-2

A method for simplifying a system of DAEs, having an index of 2 orhigher, comprises simplifying the system by projecting the system ontonormal and tangential spaces of the constraint manifold, and solving theprojected system symbolically so as to reduce the number of differentialand algebraic variables. This may be achieved using a computer, forexample a computer including a processing unit executing a Hessenbergtransforming program operable to transform the DAEs into Hessenbergform.

Systems having an index of 2 are generally expressed in Hessenberg formas:{dot over (x)}=f(t,x,z)0=h(t,x)  (2)where x, z, f, and h are defined as above for the index-1 form ofEquation 1, and the product of Jacobians

$\frac{\partial h}{\partial x}\frac{\partial f}{\partial z}$is nonsingular for all t.

Reducing system (2) to an index of 1 requires the differentiating theexpression once, so as to obtain the equation

$\begin{matrix}{0 = {{\frac{\partial{h( {t,x} )}}{\partial x}\overset{.}{x}} + \frac{\partial{h( {t,x} )}}{\partial t}}} & (3)\end{matrix}$

Accordingly, the reduced system (2) may be rewritten with equation (3)so as to obtain:

$\begin{matrix}{{\overset{.}{x} = {f( {t,x,z} )}},{0 = {{\frac{\partial{h( {t,x} )}}{\partial x}{f( {t,x,z} )}} + \frac{\partial{h( {t,x} )}}{\partial t}}}} & (4)\end{matrix}$

System (4) is in reduced form, meaning it is now in Hessenberg index-1.Simplification of system (2) begins by projecting system (4) onto thespace normal and tangential to the constraint manifold. 0=h(t,x) definesa time varying constraint manifold M(t) in the n-dimensional vectorspace of the differential variables. In order for the differentialvariables to stay on the constraint manifold for every t, {dot over (x)}must be tangent to M(t) for every t. The matrix whose rows span thenormal space of the constraint manifold is referenced herein as “C.”Since the Hessenberg form has an index of 2,

$\frac{\partial h}{\partial x}\frac{\partial f}{\partial z}$is nonsingular and the matrix “C” is of maximal row-rank. Therelationship may be expressed as

$C = {\frac{\partial h}{\partial x}.}$Differentiating the constraint manifold with respect to time yields

${{C\overset{.}{x}} + \frac{\partial h}{\partial t}} = 0.$

Describing the projection of {dot over (x)} onto the tangent space maybe done by introducing a matrix D such that D is of maximal column-rankand CD=0. The columns of D define the space which is orthogonal to thespace spanned by the rows of C. Thus, for any t, the matrix C and Dprovide two subspaces of R^(n), spanned by C and D, that are orthogonalcomplements of each other. By the theorem on orthogonal decomposition,any element R^(n) can be represented as a sum of projections onto thesesubspaces. Consequently, for any t, the derivatives of differentialvariables may be written as:{dot over (x)}=Du+C ^(T) v  (5)for certain vectors uεR^(k), wherein “u” is the tangential component and“v” is the normal component of {dot over (x)}.

Accordingly, the normal component “v” of {dot over (x)} may be obtainedby differentiating with respect to time the constraint manifold and thenprojecting the result onto the normal space spanned by the rows of C,achieving the following expression:

$\begin{matrix}{v = {{- ( {CC}^{T} )^{- 1}}\frac{\partial h}{\partial t}}} & (6)\end{matrix}$

The normal component “u” of {dot over (x)} may be obtained by projectingthe DAE ({dot over (x)}=f(t,x,z)) onto the tangential space spanned bythe columns of D and substituting the representation of {dot over (x)}as a sum of the projections onto the normal and tangential spaces of themanifold to achieve the following expression:u=(D ^(T) D)⁻¹ D ^(T) f(t,x,z)  (7)

Substituting the expressions for the normal and tangential components(5) and (6) into system (4) yields:

$\begin{matrix}{{\overset{.}{x} = {{{D( {D^{T}D} )}^{- 1}D^{T}{f( {t,x,z} )}} + {C^{T}( {{- ( {CC}^{T} )^{- 1}}\frac{\partial h}{\partial t}} )}}},{0 = {{{Cf}( {t,x,z} )} + {\frac{\partial h}{\partial t}.}}}} & (8)\end{matrix}$where the first subsystem of the system is decomposed into the sum oftangent and normal components.

Accordingly, the decomposed system (8) has the same number of equationsand variables as the originally expressed DAE. However, whenf(t,x,z)=C^(T)φ(t,x,z) for a function φ, f may be represented by aproduct of an arbitrary (n-k)-dimensional vector-valued functiondependent on t, x, and z and the transpose of the matrix C, and thesystem (8) may be simplified, meaning the number of algebraic equationsand differential variables may be reduced. Specifically, substitutingf(t,x,z)=C^(T)φ(t,x,z) reduces system (8) to:{dot over (x)}=C ^(T) v,0=φ(t,x,z)−v  (9)where Equation 9(b) can be solved with respect to z, and v is given by(6). Accordingly, projection simplifies system (2) when the aboveconditions are apparent. Specifically, system (4) has been split into asystem of ordinary differential equation and a system of algebraicequations.

Though the above case was very specific, projection may simplify abroader class of system (2). Specifically, projection of system (2)where at least one of the algebraic constraints in the vector ofdifferential variables is linear, projecting the system onto the normaland tangential spaces spanning the constraint manifold will reduce thenumber of algebraic variables and differential variables of the system,wherein the system is of the form:{dot over (x)}=f(t,x,z),0=C ₁(t)x+α(t),0=h ₂(t,x),  (10)where C₁ is an l×n matrix of maximal row-rank, meaning that l≦n, and h₂is a (k−l)-dimensional vector function which is nonlinear in x, where las used herein refers to the number of rows of C₁ and n refers to thenumber of columns of C₁.

We introduce a matrix D₁ of maximal column-rank such that C₁D₁=0. Thus,as the rows of C₁ span a subspace in R^(N) for every t, it follows thatthe columns of D₁ span the orthogonal complement to the subspace spannedby C₁, and any xεR^(N) can be represented asx=D ₁ χ+C ₁ ^(T)Ψ,  (11)where χ=(χ₁, . . . , χ_(n-1))^(T) represents projection of x onto thetangential space of the manifold formed by the linear constraints, andΨ=Ψ₁, . . . , Ψ₁)^(T) represents projection onto normal space of themanifold formed by the linear constraints. Substituting (11) into0=C₁(t)x+α(t) of (10) and making use of the relation C₁D₁=0, Ψ may besolved for, where Ψ=−(C₁C₁ ^(T))⁻¹α. From (11) we can solve for {dotover (x)}, wherein{dot over (x)}+{dot over (D)} ₁ χ+D ₁ {dot over (χ)}+Ċ ₁ ^(T) Ψ+C ₁^(T){dot over (Ψ)}  (12)

The differential equation for {dot over (χ)}, the only unknown part of{dot over (x)} after projection, may be derived by substituting (12)into {dot over (x)}=f(t,x,z) and left-multiplication by D₁ ^(T), where{dot over (χ)}=(D ₁ ^(T) D ₁)⁻¹ D ₁ ^(T) [f(t,D ₁ χ+C ₁ ^(T) Ψ,z)−{dotover (D)} ₁ χ−Ċ ₁ ^(T)Ψ]  (13)

In expression (13), z is unknown. However, z may be found bydifferentiating algebraic equations 0=C₁(t)x+φ(t), and 0=h₂(t,x), so asto obtain

$\begin{matrix}{{{{{C\overset{.}{x}} + \eta} = 0},{where}}{C = {\begin{pmatrix}C_{1} \\\frac{\partial h_{2}}{\partial x}\end{pmatrix}\mspace{14mu}{and}}}{\eta = {\begin{pmatrix}{{{\overset{.}{C}}_{1}x} + \overset{.}{a}} \\\frac{\partial h_{2}}{\partial t}\end{pmatrix}.}}} & (14)\end{matrix}$On the other hand, left-multiplying {dot over (x)}=f(t,x,z) by C yields:C{dot over (x)}=Cf  (15)

Using equations (14) and (5), an expression for z may be obtained wherez is solvable because

$\begin{pmatrix}C_{1} \\\frac{\partial h_{2}}{\partial x}\end{pmatrix}\frac{\partial f}{\partial z}$is a nonsingular matrix. Specifically, the combination of equations (13)and (14) yields the expression:Cf(t,D ₁ χ+C ₁ ^(T) Ψ,z)+η(t,D ₁ χ+C ₁ ^(T)Ψ)=0  (16)

Thus the DAE having n number of differential variables and k number ofalgebraic variables is simplified to the DAE system described inequations (13) and (16), which has n−l differential variables, where lis the number of linear constraints, and k is a number of algebraicvariables.

Accordingly, projection simplifies a system, where the system is in theform (10), and where Ψ=−(C₁C₁ ^(T))⁻¹α. Using the relationship C₁D₁=0,D₁ may be found symbolically, and thus simplified DAEs (13) and (16) maybe numerically solved for χ and z, allowing x to be known using Ψ and χin x=D₁χ+C₁ ^(T)Ψ.

Hessenberg Index-3

Simplification of a DAE having a higher index is also possible usingprojection. It is known that a Hessenberg index-3 may be generallyexpressed as follows:{dot over (y)}=f(t,x,y,z),{dot over (x)}=g(t,x,y),0=h(t,x)  (17)where y=y(t)=(y₁(t), . . . y_(m)(f))^(T) and x=x(t)=(x₁(t), . . .x_(n)(t))^(T) are the vectors of differential variables, f=(f₁, . . . ,f_(m))^(T), g=(g, . . . , g_(n))^(T) are known vector functions,representing differential equations, z=z(t)=z₁(t), . . . z_(k)(t)^(T) isthe vector of algebraic variables, h=(h₁, . . . , h_(k))^(T) is a knownvector function representing algebraic equations and the product of theJacobians

$\frac{\partial h}{\partial x}\frac{\partial g}{\partial y}\frac{\partial f}{\partial z}$is nonsingular for all t.

Thus, {dot over (y)}=f(t,x,y,z) (Equation 17a) and {dot over(x)}=g(t,x,y) (17b) are two different types of differential equations,and 0=h(t,x)(17c) is a set of algebraic equations defining theconstraint manifold. Reduction of the Hessenberg index-3 form to index 1consists of two differentiations of the system (17) and projecting thesecond differentiation onto the normal and tangential spaces of theconstraint manifold.

Specifically, {dot over (x)}=g(t,x,y), and 0=h(t,x) of (17) aredifferentiated and then projected. As before, projection is achieved byintroducing the tangential and normal spaces of the constrained manifoldso as to yield expressions (5) and (6)

Differentiating {dot over (x)}=g(t,x,y), yields

$\begin{matrix}{\overset{¨}{x} = {\frac{\partial g}{\partial t} + {\frac{\partial g}{\partial x}g} + {\frac{\partial g}{\partial y}f}}} & (18)\end{matrix}$

Accordingly, using (5), (18) may be expressed as:

$\begin{matrix}{{{\overset{.}{D}u} + {D\overset{.}{u}} + {{\overset{.}{C}}^{T}v} + {C^{T}\overset{.}{v}}} = {\frac{\partial g}{\partial t} + {\frac{\partial g}{\partial x}g} + {\frac{\partial g}{\partial y}f}}} & (19)\end{matrix}$

Next, expression (19) is also projected onto the tangential space. Thismay be done by left-multiplying expression (19) by D^(T), resulting in

$\begin{matrix}{{{D^{T}\overset{.}{D}u} + {D^{T}D\overset{.}{u}} + {D^{T}{\overset{.}{C}}^{T}v}} = {{D^{T}\frac{\partial g}{\partial t}} + {D^{T}\frac{\partial g}{\partial x}g} + {D^{T}\frac{\partial g}{\partial y}f}}} & (20)\end{matrix}$

Since D is of maximal column-rank, we can express {dot over (u)} interms of t, x, y, and u. Thus, “v” is known from (6), and u may bederived from differential equation (20).

However, the term “y” remains unknown, which may be referred to as alevel-1 differential variable, described further below, and may beeliminated if represented in terms of t, x, and u. This may be done bysubstituting the representation (5) of {dot over (x)} into (17b) {dotover (x)}=g(t,x,y), accordingly it follows that:{dot over (x)}=Du+C ^(T) v=g(t,x,y)  (21)

Assuming that the Jacobian ∂g/∂y is invertible for any t, and the numberof components in vector x is equal to the number of components in y,expression (21) may be solved uniquely for y. Thus, the only unknownvariable remaining is z. Later, a generalized approach is introducedthat does not require these assumptions.

In order to reduce the Hessenberg form to an index 1, the methodproceeds to the step of obtaining an algebraic equation solvable for z.The algebraic equation may be obtained by projecting (19) onto thenormal space. In order to obtain Ċ and {dot over (v)}, system (17c) isdifferentiated a second time. Thus, the following index 1 DAEs may beobtained utilizing the expressions (6), (19), (20), and (21), asfollows:

$\begin{matrix}{{\overset{.}{x} = {{Du} + {C^{T}v}}},{\overset{.}{u} = {( {D^{T}D} )^{- 1}{D^{T}\lbrack {\frac{\partial g}{\partial t} + {\frac{\partial g}{\partial x}g} + {\frac{\partial g}{\partial y}f} - {\overset{.}{D}u} - {{\overset{.}{C}}^{T}v}} \rbrack}}},{0 = {{Du} + {C^{T}v} - g}},{0 = {{\overset{.}{C}g} + {C\lbrack {\frac{\partial g}{\partial t} + {\frac{\partial g}{\partial x}g} + {\frac{\partial g}{\partial y}f}} \rbrack} + \frac{\partial^{2}h}{\partial t^{2}} + {\frac{\partial^{2}h}{{\partial x}{\partial t}}g}}}} & (22)\end{matrix}$where “v” is given by expression (6). The projected expression of thereduced Hessenberg form has 2n−k differential variables, “x” and “u”,and n+k algebraic variables, “y” and “z”.

As with the case above wherein projection may be used to reduce thenumber of differential and algebraic variables in the system, the numberof algebraic variables and differential variables may also be reduced byprojecting system onto the normal and tangential subspaces of theconstraint manifold of the system where at least one algebraicconstraint in the vector of differential variables is linear.

A method pertaining to a DAE having an index of 3, when some constraintsare linear, is now described where the Hessenberg index-3 system is ofthe form:{dot over (y)}=f(t,x,y,z),{dot over (x)}=g(t,x,y),0=C ₁(t)x+α(t),0=h ₂(t,x),  (23)where C₁ is an l×n matrix of maximal row rank, such that l≦n, and h₂ isa (k−l)-dimensional vector function which is nonlinear in x.

Subsystems {dot over (x)}=g(t,x,y), 0=C₁(t)x+α(t), and 0=h₂(t,x) may berewritten to take into consideration the normal and tangential componentof the vector of differential variables, as expressed in equation (11),where the vector Ψ is given by Ψ=−(C₁C₁ ^(T))⁻α. The equation for χ,derived from (13), is:{dot over (χ)}=(D ₁ ^(T) D ₁)⁻¹ D ₁ ^(T) [g(t,D ₁ χ+C ₁ ^(T) Ψ,y)−{dotover (D)} ₁ χ−Ċ ₁ ^(T)Ψ]  (24)

Accordingly, it follows that subsystem (22a) {dot over (χ)} may bereplaced with (24). Thus system (23) may be reduced to:

$\begin{matrix}{{\overset{.}{\chi} = {( {D_{1}^{T}D_{1}} )^{- 1}{D_{1}^{T}\lbrack {g - {\overset{.}{D}}_{1\chi} - {{\overset{.}{C}}_{1}^{T}\psi}} \rbrack}}},{\overset{.}{u} = {( {D^{T}D} )^{- 1}{D^{T}\lbrack {\frac{\partial g}{\partial t} + {\frac{\partial g}{\partial x}g} + {\frac{\partial g}{\partial y}f} - {\overset{.}{D}u} - {{\overset{.}{C}}^{T}v}} \rbrack}}},{0 = {{Du} + {C^{T}v} - g}},{0 = {{\overset{.}{C}g} + {C\lbrack {\frac{\partial g}{\partial t} + {\frac{\partial g}{\partial x}g} + {\frac{\partial g}{\partial y}f}} \rbrack} + \frac{\partial^{2}h}{\partial t^{2}} + {\frac{\partial^{2}h}{{\partial x}{\partial t}}g}}},} & (25)\end{matrix}$Accordingly, system (25) has (2n−k−l)-differential variables χ and u andn+k algebraic variables y and z.

Projection may also simplify a Hessenberg Index-3 DAE where thesubsystem, 0=Du+C^(T)v−g in system (22) or (25) is solvable for “y”symbolically.

In the Hessenberg Index-3 case, there are two different classes ofdifferential variables, represented by variables “x” and “y.” Thevariables “x” represent differential variables that appear in thealgebraic equations of the system and such that there are no algebraicvariables in the differential equations for “x.” The variables “x”represent the variables that are affected by the constraints. Thederivatives of the other class of differential variables, “y”, representthe means by which the constraint are enforced.

For example, for a mechanical system subject to holonomic constraints,the constraints are expressed in coordinates, variables x, and areenforced by the constraint forces. By the Newton's law, the forces areconnected to the acceleration, which plays the role of {dot over (y)}.Thus, y represents the velocities. If Equation 22c or 25c may be solvedsymbolically for “y” then “y” may be removed from system (22) or (25),respectively, so as to decouple the system. Accordingly, the remainingsubsystems of the system form a DAE system with respect to the variables“x”, “u”, and “z”. Thus, n algebraic equations are removed from thesystem.

Additionally, projection may also simplify Hessenberg Index-3 DAEs wherethe subsystem,

${0 = {{\overset{.}{C}g} + {C\lbrack {\frac{\partial g}{\partial t} + {\frac{\partial g}{\partial x}g} + {\frac{\partial g}{\partial y}f}} \rbrack} + \frac{\partial^{2}h}{\partial t^{2}} + {\frac{\partial^{2}h}{{\partial x}{\partial t}}g}}},$(22d) or (25d), is solvable symbolically for “z.”

If subsystem (22d) or (25d) is solvable for “z” symbolically, thensubsystems (22a-22c) or (25a-25c) become DAE systems with respect to“x”, “u” and “y” (or “x”, “u” and “y”), respectively. By solving thesubsystem (22d or 25d) symbolically and substituting the result into(22a-22c) or (25a-25c), respectively, k algebraic equations and kalgebraic variables may be removed from the corresponding system.

The Projection Method and Higher Hessenberg Index DAEs

We generalize the approach from the previous section to DAEs ofarbitrary Hessenberg index. We consider a Hessenberg index-(r+2) DAE

$\begin{matrix}{{{\overset{.}{y}}_{r} = {f_{r}( {t,x,y_{1},y_{2},\ldots\mspace{14mu},y_{r},z} )}},\vdots} & ( {26a} ) \\{{{\overset{.}{y}}_{2} = {f_{2}( {t,x,y_{1},y_{2},y_{3}} )}},} & ( {26b} ) \\{{{\overset{.}{y}}_{1} = {f_{1}( {t,x,y_{1},y_{2}} )}},} & ( {26c} ) \\{{\overset{.}{x} = {g( {t,x,y_{1}} )}},} & ( {26d} ) \\{{0 = {h( {t,x} )}},{where}} & ( {26e} ) \\{\frac{\partial h}{\partial x}\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}}\mspace{14mu}\ldots\mspace{14mu}\frac{\partial f_{r - 1}}{\partial y_{r}}\frac{\partial f_{r}}{\partial z}{\mspace{11mu}\;}{is}\mspace{14mu}{nonsingular}\mspace{14mu}{for}\mspace{14mu}{all}\mspace{14mu}{t.}} & (27)\end{matrix}$

The vectors y₁, y₂, . . . , y_(r) are of size m₁, m₂, . . . , m_(r),respectively. Due to condition (27), the ranks of the Jacobins

$\frac{\partial g}{\partial y_{1}},\frac{\partial f_{1}}{\partial y_{2}}\;,\ldots\mspace{14mu},\;\frac{\partial f_{r - 1}}{\partial y_{r}}$and their products forming subproducts in condition (27) are at least kat any t. We re-order the components of the vectors y₁, y₂, . . . ,y_(r) so that the first j₁, j₂, . . . , j_(r) components, respectively,denoted by superscript a, are such that

$\begin{matrix}{{{{rank}\frac{\partial g}{\partial y_{1}}} = {{{rank}\frac{\partial g}{\partial y_{1}^{a}}} = j_{1}}},} & ( {28a} ) \\{{{{rank}\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}}} = {{{rank}\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}^{n}}} = j_{2}}},\vdots} & ( {28b} ) \\{{{rank}\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}}\mspace{14mu}\ldots\mspace{14mu}\frac{\partial f_{r - 1}}{\partial y_{r}}} = {{{rank}\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}}\mspace{14mu}\ldots\mspace{14mu}\frac{\partial f_{r - 1}}{\partial y_{r}^{n}}} = {j_{r}.}}} & ( {28c} )\end{matrix}$

As in the previous cases the projection method leads to therepresentation

$\begin{matrix}{{\overset{.}{x} = {{Du}_{1} + {C^{T}v}}},{{{where}\mspace{14mu} C} = \frac{\partial h}{\partial x}}} & (29)\end{matrix}$is of maximal row-rank, D is of maximal column-rank and satisfies therelation

${{CD} = 0},{v = {{- ( {CC}^{T} )^{- 1}}\frac{\partial h}{\partial t}}},$and u₁ is an (n−k)-dimensional vector. Then we introduce(n−k)-dimensional vectors u₂, . . . , u_(r) such that

$\begin{matrix}{{{\overset{.}{u}}_{1} = u_{2}},} & ( {30a} ) \\{{{\overset{.}{u}}_{2} = u_{3}},\vdots} & ( {30b} ) \\{{{\overset{.}{u}}_{r - 1} = u_{r}},} & ( {30c} )\end{matrix}$and we want to re-write equations (26) in terms of the projection of{dot over (x)} onto the tangential space, u₁, and its derivatives. Wesubstitute equation (29) into equation (26d), obtainingDu ₁ +C ^(T) v=g(t,x,y ₁).  (31)From condition (28a), it follows that equation (31) can be solved for y₁^(a).

Differentiating equations (31) with respect to time, we obtain

$\begin{matrix}{{{{\overset{.}{D}u_{1}} + {Du}_{2} + {{\overset{.}{C}}^{T}v} + {C^{T}\overset{.}{v}}} = {\frac{\partial g}{\partial t} + {\frac{\partial g}{\partial x}g} + {\frac{\partial g}{\partial y_{1}}f_{1}}}},} & (32)\end{matrix}$where we made use of equations (26d), (26c), and (30a). From condition(28b), it follows that equation (32) can be solved for y₂ ^(a).

Differentiating equation (31) with respect to time second time, i.e.differentiating equation (32), and replacing derivatives of variablesthat occur in the left-hand side with their representation according toequations (29) and (30) and derivatives that occur in the right-handside according to equations (26), we obtain algebraic equations that canbe solved for y₃ ^(a). We continue this process until we differentiateequation (31) r times. After each differentiation we make use ofequations (29) and (30) in the left-hand side and of equations (26) onthe right-hand side. The i-th differentiation gives us an algebraicequation, which due to conditions (28) can be solved for y_(i+1) ^(a),for i=1, . . . , r−1.

The r-th differentiation gives us an equation, from which we derive adifferential equation for {dot over (u)}_(r) and an algebraic equationfor z. Let us consider the result of r-th differentiation. We areinterested only in the terms that define equations for {dot over(u)}_(r) and z, and, therefore, for simplicity, we represent the resultof r-th differentiation as follows:

$\begin{matrix}{{{{D{\overset{.}{u}}_{r}} + {\varphi_{1}( {t,x,u_{1},u_{2},\ldots\mspace{14mu},u_{r}} )}} = {{\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}}\frac{\partial f_{2}}{\partial y_{3}}\mspace{14mu}\ldots\mspace{20mu}\frac{\partial f_{r - 1}}{\partial y_{r}}f_{r}} + {\varphi_{2}( {t,x,y_{1},y_{2},\ldots\mspace{14mu},y_{r}} )}}},} & (33)\end{matrix}$where functions φ₁ and φ₂ represent the remaining terms, in which we arenot interested.

If we project equation (33) onto the tangential space of the constraintmanifold by multiplying both sides of the equation by D^(T) on the left,we obtain an equation for {dot over (u)}_(r).

$\begin{matrix}{{\overset{.}{u}}_{r} = {( {D^{T}D} )^{- 1}{{D^{T}\lbrack {{\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}}\frac{\partial f_{2}}{\partial y_{3}}\mspace{14mu}\ldots\mspace{20mu}\frac{\partial f_{r - 1}}{\partial y_{r}}f_{r}} + \varphi_{2} - \varphi_{1}} \rbrack}.}}} & (34)\end{matrix}$

if we project equation (33) onto the normal space of the constraintmanifold by multiplying both sides of the equation by C on the left, weobtain an equation for z, which is solvable due to condition (27) and isindependent of {dot over (u)}_(r) due to relation CD=0.

$\begin{matrix}{0 = {( {CC}^{T} )^{- 1}{{C\lbrack {{\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}}\frac{\partial f_{2}}{\partial y_{3}}\mspace{14mu}\ldots\mspace{20mu}\frac{\partial f_{r - 1}}{\partial y_{r}}f_{r}} + \varphi_{2} - \varphi_{1}} \rbrack}.}}} & (35)\end{matrix}$

The DAE system simplified by the projection method is then given byequations (29), (30), (31), (32), (34), (35), algebraic equations for y₃^(a), . . . , y_(r) ^(a) and differential equations for y₁ ^(d), . . . ,y_(r) ^(d), where we represent each y_(i) in terms of y_(i) ^(a) andy_(i) ^(d), for i=1, . . . , r.

$\begin{matrix}{{\overset{.}{x} = {{Du}_{1} + {C^{T}v}}},} & ( {36a} ) \\{{{\overset{.}{u}}_{1} = u_{2}},} & ( {36b} ) \\{{{\overset{.}{u}}_{2} = u_{3}},\vdots} & ( {36c} ) \\{{{\overset{.}{u}}_{r - 1} = u_{r}},} & ( {36d} ) \\{{{\overset{.}{u}}_{r} = {( {D^{T}D} )^{- 1}{D^{T}\lbrack {{\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}}\frac{\partial f_{2}}{\partial y_{3}}\mspace{14mu}\ldots\mspace{20mu}\frac{\partial f_{r - 1}}{\partial y_{r}}f_{r}} + \varphi_{2} - \varphi_{1}} \rbrack}}},} & ( {36e} ) \\{{{\overset{.}{y}}_{1}^{d} = f_{1}^{d}},} & ( {36f} ) \\{{{\overset{.}{y}}_{2}^{d} = f_{2}^{d}},\vdots} & ( {36g} ) \\{{{\overset{.}{y}}_{r}^{d} = f_{r}^{d}},} & ( {36h} ) \\{{0 = {{Du}_{1} + {C^{T}v} - g}},} & ( {36i} ) \\{{0 = {{\overset{.}{D}u_{1}} + {Du}_{2} + {{\overset{.}{C}}^{T}v} + {C^{T}\overset{.}{v}} - \frac{\partial g}{\partial t} - {\frac{\partial g}{\partial x}g} - {\frac{\partial g}{\partial y_{1}}f_{1}}}},} & ( {36j} ) \\{{0 = {\frac{\mathbb{d}^{3}}{\mathbb{d}t^{3}}( {{Du}_{1} + {C^{T}v}} ){_{{(29)},{(30)}}{{- \frac{\mathbb{d}^{3}}{\mathbb{d}t^{3}}}(g)}}_{(26)}}},\vdots} & ( {36k} ) \\{{0 = {\frac{\mathbb{d}^{r - 1}}{\mathbb{d}t^{r - 1}}( {{Du}_{1} + {C^{T}v}} ){_{{(29)},{(30)}}{{- \frac{\mathbb{d}^{r - 1}}{\mathbb{d}t^{r - 1}}}(g)}}_{(26)}}},} & ( {36l} ) \\{0 = {( {CC}^{T} )^{- 1}{{C\lbrack {{\frac{\partial g}{\partial y_{1}}\frac{\partial f_{1}}{\partial y_{2}}\frac{\partial f_{2}}{\partial y_{3}}\mspace{14mu}\ldots\mspace{20mu}\frac{\partial f_{r - 1}}{\partial y_{r}}f_{r}} + \varphi_{2} - \varphi_{1}} \rbrack}.}}} & ( {36m} )\end{matrix}$

System (36) consists of differential equations (36a)-(36h) with respectto n+r(n−k)+(m₁−j₁)+(m₂−j₂)+ . . . +(m_(r)−j_(r)) variables x, u₁, . . ., u_(r), y₁ ^(d), . . . , y_(r) ^(d) and of algebraic equations(36i)-(36m) with respect to j₁+j₂+ . . . +j_(r)+k variables y₁ ^(a), . .. , y_(r) ^(a), z. Note that index reduction alone yields an index-1 DAEwith respect to n+m₁+m₂+ . . . +m_(r) differential and k algebraicvariables.

FIGS. 1A-1B illustrate a general approach for simplifying a DAE of indexof at least 3, in flow chart form. The DAE is received as an input(block 10), and converted to Hessenberg form and separated into threegroups (block 12), including algebraic equations, level-zero (level-0)differential equations, and differential equations of levels 1 throughr. Block 14 corresponds to separating the differential equations oflevels 1 through r. In this example, level-1 corresponds to Equation 26cand level-r corresponds to Equation 26a. Block 16 corresponds toseparating the differential equations of level 0, corresponding toEquation 26d, and block 18 corresponds to separating the algebraicmanifold equations corresponding to Equation 26e.

The manifold equations, separated at block 18, are differentiated (block20) and used to obtain an equation for the normal projection of x′ ontoa manifold v (block 22), for which v and bases for the normal andtangent spaces are determined (block 24), allowing determination of v asa function of t (block 26). This allows x′ to be represented in terms ofthe projections onto normal (v) and tangent (u) spaces, indicated atblock 28 (Equation 29). Formulae, denoted *, are then determined for x′in terms of unknown u and known v (block 30), the latter being aprojection onto normal space that is used in the simplified DAE (block64).

The equations for x′ in terms of unknown u are substituted into thedifferential equations of level zero (block 32), giving an equation ofthe form shown in block 34 (Equation 31). The value of i is assigned to1 (block 36), i is compared to r (block 38), and if i≦r, the formulaeA^(i) (block 40) are differentiated. Formulae from block 30 andequations from block 14 are substituted into the result of thedifferentiation of formulae A^(i), giving formulae A^(i+1). The value ofi is incremented by 1 (block 48), and compared to r again (block 38).

If i>r, the formulae A^(r+1) (Equation 33) obtained (block 42) areprojected onto tangent space (block 58), giving equations for the r-thderivative of u (Equations 36e, the same as equations in block 62), andused in the simplified DAE (block 64). The equations A^(r+1) obtained(block 42) are also projected onto normal space (block 56), givingalgebraic equations for z (block 60, Equation (36m)), also used in thesimplified DAE (block 64).

The components of y_(i) that can be found from equations A^(i) arecalled components y_(i) ^(a), and equations A^(i) are algebraicequations for y_(i) ^(a). The remaining components of y_(i) are calledy_(i) ^(d) and we use the corresponding original differential equationsfor these components. In this way, the formulae A^(i) (block 40) forwhich components of y_(i) (A^(i)) can be solved algebraically (block 50)are identified, giving differential equations for y_(i) ^(d) (block 52,Equations 36f-36h) and algebraic equations for y_(i) ^(a) (block 54,Equations 36i-36l), and these are further used in the formulation of thesimplified DAE (block 64).

Hence, a general approach includes transforming the differentialalgebraic equations into Hessenberg form, the Hessenberg form includinga set of algebraic equations representing constraint manifold, and a setof differential equations of levels 0 through r, where r is greater orequal to 1. Algebraic equations of the manifold are differentiated so asto produce bases for normal and tangent spaces of the manifold and anequation for projection of the derivatives of 0th level differentialvariables onto the normal space. The derivatives {dot over (x)} of the0th level differential variables x are represented in terms of theprojections onto the normal and tangent spaces (e.g. see Equation 29),substituting this representation whenever derivatives of 0th leveldifferential variables are encountered so as to convert the 0th leveldifferential equations into the 1st level algebraic equations (e.g.equations for A¹ in FIG. 1, block 34, and see also Equation 31 for the1st level differential equations).

Level-1 differential variables y₁ for which the level-1 algebraicequations can and cannot be solved are identified, and only the level-1differential equations for those level 1 differential variables forwhich level-1 algebraic equations cannot be solved are retained, so asto reduce the number of differential variables.

Level-1 algebraic equations are differentiated r−1 times, each timesubstituting differential equations of levels 1 through r−1, so as toobtain algebraic equations of levels 2 through r (e.g. equations (A^(i))in the loop of the FIG. 1 for i=1 . . . r, the loop including blocks 38,40, 44, 46, and 48) with respect to differential variables of levels 2through r, respectively.

Differential variables of each level from 1 to r−1, y_(t), . . .y_(r−1), for which algebraic equations of the corresponding level canand cannot be solved are identified, and differential equations of thecorresponding level for those variables for which algebraic equations ofthe corresponding level cannot be solved are kept, so as to reduce thenumber of differential variables.

Level-1 algebraic equations are then differentiated the r-th time (e.g.yielding Equation 33), substituting differential equations of levels 1through r, and projecting the result onto normal and tangent spaces soas to obtain level-(r+1) algebraic equations for algebraic variables(e.g. block 60 in FIG. 1), and equations for the r-th derivative of thetangent projection of the derivative of level-0 differential variables,respectively, and thus to reduce the index of the DAE and number ofdifferential variables.

For a method such as described above, if the level-0 algebraic equationsare linear with respect to at least one of the level-0 differentialvariables, the method may further include considering a linear manifolddefined by linear algebraic equations only, representing level-0differential variables in terms of the projection to the linear manifold(its tangent space is the manifold itself) and to its normal space,using linear algebraic equations to find the normal projection of thelevel-0 differential variables, and rewriting differential equations oflevel 0 as differential equations for the projection of the level-0differential variables onto the manifold so as to reduce the number ofdifferential variables.

If algebraic equations of level i can be solved symbolically for some ofthe differential variables of level i, the method may further include astep of solving these equations and substituting them into the system soas to reduce the number of algebraic variables in the system.

If the algebraic equations of level r+1 can be solved symbolically forsome of the algebraic variables, the general method may further includea step of solving these equations and substituting them into the systemso as to reduce the number of algebraic variables in the system.

An example method of simplifying a system expressed as differentialalgebraic equations includes providing representations of the algebraicequations as a constrained manifold and determining the normal andtangential spaces of the manifold, transforming the DAE into Hessenbergform, reducing the Hessenberg form to an index of 1, and projecting thesystem onto the tangential and normal spaces of the manifold so as toreduce the number of algebraic and differential variables in the system.

Examples of the present invention may be implemented in a virtualengineering modeling environment, which may include a computer systemcomprising one or more of the following: a processor, memory component,user interface, input/output device, data buses providing communicationsbetween various components, and the like. A modeling environment mayinclude a network of computers. An initial configuration of the modelingenvironment, such as a computer system, may include or receive aninitial model of a physical system in terms of a DAE. The computersystem may further allow test and verification of the physical systemwithin the virtual engineering environment.

Examples of the present invention allow improved system modeling andvirtual prototyping by replacing the initial model with a revised modelin terms of ODEs, that may include fewer differential and/or algebraicparameters. This allows significantly faster computation of modelpredictions and behavior, and hence reduced prototyping time and fasterproduct development.

Examples of the present invention include methods performed by acomputer algebra system capable of symbolic mathematics, such as avirtual engineering environment. A virtual engineering environment maycomprise one or more processors, memory components such as volatileand/or non-volatile memory components, and one or more communicationsinterfaces such as a display, data entry components, and the like. Avirtual engineering computer may be operable to numerically solve a DAEmodel of a physical system using approaches described herein.

The invention is not restricted to the illustrative examples describedabove. Examples described are not intended to limit the scope of theinvention. Changes therein, other combinations of elements, and otheruses will occur to those skilled in the art.

Modifications and variations of the present invention are possible inlight of the above teachings and may be practiced otherwise than asspecifically described while within the scope of the appended claims.

We claim:
 1. In a computational environment having at least onecomputing system including at least one processor, a method ofprototyping and testing a virtual model representing a physical system,wherein the method includes simplifying the model defined by a system ofdifferential algebraic equations (DAEs) and not being in Hessenbergform, the method comprising: transforming on at least one computingsystem the DAEs into Hessenberg form defining an index and having levelsfrom level-0 through level-r, wherein r is a variable representing alevel greater than zero, the Hessenberg form including algebraicequations representing a constraint manifold, and differentialequations; differentiating on at least one computing system thealgebraic equations to produce bases for a normal space and a tangentspace of the constraint manifold; representing derivatives on at leastone computing system of level-zero differential variables as projectionsonto the normal space and the tangent space of the constraint manifold,and substituting the projections into the level-zero differentialequations, so as to reduce the number of differential variables andalgebraic variables in the model of the physical system; converting onat least one computing system the level-zero differential equations intolevel-1 algebraic equations by representing derivatives of thelevel-zero differential variables in terms of the projections onto thenormal and tangent spaces; retaining on at least one computing systemonly the level-1 differential equations corresponding to level-1differential variables for which level-1 algebraic equations cannot besolved; differentiating on at least one computing system the level-1algebraic equations r−1 times, each time substituting differentialequations of levels 1 through r−1, so as to obtain algebraic equationsof levels 2 through r with respect to differential variables of levels 2through r, respectively, keeping only differential equations for whichalgebraic equations of the corresponding level cannot be solved;differentiating on at least one computing system the level-1 algebraicequations again, substituting differential equations of levels 1 throughr to generate a result, and projecting the result onto normal andtangent spaces so as to obtain level-(r+1) algebraic equations foralgebraic variables and equations for the r-th derivative of the tangentprojection of the derivative of level-0 differential variables,respectively; wherein the model of the physical system is simplified byreducing the index of the DAE and by reducing the number of differentialvariables in the model of the physical system; and wherein thesimplified model of the physical system generates parameters fordesigning of a prototype and virtual testing effecting the physicalsystem determined by variation of variables of the model.
 2. The methodof claim 1, wherein level-zero algebraic equations are linear withrespect to at least one of the level-zero differential variables, themethod further including: establishing a linear manifold defined bylinear algebraic equations only; representing level-zero differentialvariables in terms of projections to the linear manifold and to normalspace; using linear algebraic equations to find the normal projection ofthe level-zero differential variables; and rewriting level-zerodifferential equations as differential equations for the projection ofthe level-zero differential variables onto the linear manifold, so as toreduce the number of differential variables in the model of the physicalsystem.
 3. The method of claim 1, wherein the algebraic equations of adefined level are solved symbolically for at least one of thedifferential variables of the defined level, the method furtherincluding solving the algebraic equations so as to reduce the number ofalgebraic variables in the system.
 4. The method of claim 3, whereinalgebraic equations of level r+1 are solved symbolically for at leastone of the algebraic variables, the method further including solving thealgebraic equations so as to reduce the number of algebraic variables inthe system.
 5. The method of claim 4, the method removing all algebraicvariables and at least one differential variable from the model of thephysical system.
 6. The method of claim 1, the DAEs having levels fromlevel-0 through level-r, where r>1.
 7. The method of claim 1, thecomputational environment being an engineering computer configured toperform symbolic mathematics.
 8. The method of claim 1, wherein thephysical system is a mechanical system.
 9. The method of claim 8 whereinthe mechanical system is an engine model.
 10. The method of claim 1wherein the physical system is an electrical system.
 11. The method ofclaim 10 wherein the electrical system is an electrical circuit model.